library(readstata13)
library(tidyverse)
library(fixest)
library(readxl)
library(broom)
library(ggthemes)
library(modelsummary)
library(ROCit)
library(pROC)
library(dplyr)

if(!require('data.table')) {
  install.packages('data.table')
  library('data.table')
}
directory = c("/Users/roniyadlin/Dropbox/Climate conflict")
nd_conflict <- read.dta13("/Users/roniyadlin/Dropbox/Climate conflict/data/clean/climate_conflict_panel.dta")

##Run original OLS regression to get coefficients
origdf<-nd_conflict %>%
  filter(year>2000)   %>%
  filter(year<2010) 

setFixest_dict(c(anyconflict_1y = "1 year", anyconflict_3y = "3 years",
                 anyconflict_5y = "5 years", anynatdis = "Any natural disaster",
                 country = "Country", year="Year"))

##One year regression
conflict_est_1yr <- feols(anyconflict_1y
                      ~ anynatdis | country[year] + year, 
                      data = origdf)
beta_1yr<-unname(conflict_est_1yr$coefficients)
fixef_1yr<-fixef(conflict_est_1yr)
int_1yr<-.0742

##Three year regression
conflict_est_3yr <- feols(anyconflict_3y
                          ~ anynatdis | country[year] + year, 
                          data = origdf)
beta_3yr<-unname(conflict_est_3yr$coefficients)
fixef_3yr<-fixef(conflict_est_3yr)
int_3yr<-.1103

#Five year regression
conflict_est_5yr <- feols(anyconflict_5y
                          ~ anynatdis | country[year] + year, 
                          data = origdf)
beta_5yr<-unname(conflict_est_5yr$coefficients)
fixef_5yr<-fixef(conflict_est_5yr)
int_5yr<-.1314

##Prediction using 2010-current
predictdf<-nd_conflict %>%
  filter(year>2009)   %>%
  filter(year<2022) 

pred_df<-data.frame(predictdf$year,predictdf$country,predictdf$anynatdis)
colnames(pred_df)<-c("Year","Country","Any Natural Disaster")
totobs<-(nrow(pred_df))

#Was there a war?
totevents<-data.frame(predictdf$total_events)
totevents$War = ifelse(totevents > 0, 1, 0)
totwars<-sum(na.omit(totevents$War))

##Prediction using 1 year regression
#Extract fixed effects from regression results
fevec_cy1yr<-data.frame(fixef_1yr$`country[[year]]`)
fevec_cy1yr <- tibble::rownames_to_column(fevec_cy1yr, "Country")


#Merge fixed effects with new data
pred_df_1yr<-merge.data.frame(pred_df,fevec_cy1yr,by="Country")

#Predict using regression equation--Does not have year fixed effects because 
#they are in the future...
pred_1yr<-int_1yr+beta_1yr%*%predictdf$anynatdis+pred_df_1yr$fixef_1yr..country..year...
pred_df_1yr<-data.frame(pred_df_1yr,t(pred_1yr))
colnames(pred_df_1yr)<-c("Country","Year","Any Natural Disaster",
                         "Country Year Fixed Effects","Prob of Conflict")

##Prediction using 3 year regression results
fevec_cy3yr<-data.frame(fixef_3yr$`country[[year]]`)
fevec_cy3yr <- tibble::rownames_to_column(fevec_cy3yr, "Country")
pred_df_3yr<-merge.data.frame(pred_df,fevec_cy3yr,by="Country")
pred_3yr<-int_3yr+beta_3yr%*%predictdf$anynatdis+pred_df_3yr$fixef_3yr..country..year...
pred_df_3yr<-data.frame(pred_df_3yr,t(pred_3yr))
colnames(pred_df_3yr)<-c("Country","Year","Any Natural Disaster",
                         "Country Year Fixed Effects","Prob of Conflict")

##Prediction using 5 year regression results
fevec_cy5yr<-data.frame(fixef_5yr$`country[[year]]`)
fevec_cy5yr <- tibble::rownames_to_column(fevec_cy5yr, "Country")
pred_df_5yr<-merge.data.frame(pred_df,fevec_cy5yr,by="Country")

pred_5yr<-int_5yr+beta_5yr%*%predictdf$anynatdis+pred_df_5yr$fixef_5yr..country..year...
pred_df_5yr<-data.frame(pred_df_5yr,t(pred_5yr))
colnames(pred_df_5yr)<-c("Country","Year","Any Natural Disaster",
                         "Country Year Fixed Effects","Prob of Conflict")
#ROC Curve
roc1<-roc(totevents$War,pred_df_1yr$`Prob of Conflict`)
gl <- ggroc(roc1, color="red", legacy.axes = TRUE)
glsave<-gl + xlab("FPR") + ylab("TPR") +
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="blue", linetype="dashed")
ggsave("/Users/roniyadlin/Dropbox/Climate conflict/output/figures/ROC_LY_1yr.png",glsave,device="png")

roc3<-roc(totevents$War,pred_df_3yr$`Prob of Conflict`)
gl <- ggroc(roc3, color="red", legacy.axes = TRUE)
glsave<-gl + xlab("FPR") + ylab("TPR") +
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="blue", linetype="dashed")
ggsave("/Users/roniyadlin/Dropbox/Climate conflict/output/figures/ROC_LY_3yr.png",glsave,device="png")

roc5<-roc(totevents$War,pred_df_5yr$`Prob of Conflict`)
gl <- ggroc(roc5, color="red", legacy.axes = TRUE)
glsave<-gl + xlab("FPR") + ylab("TPR") +
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="blue", linetype="dashed")
ggsave("/Users/roniyadlin/Dropbox/Climate conflict/output/figures/ROC_LY_5yr.png",glsave,device="png")

#Does the prediction match actual wars
pred_df_1yr$War = ifelse(pred_df_1yr$`Prob of Conflict` > .172, 1, 0)
compdf_1y<-data.frame(pred_df_1yr,totevents$War)
colnames(compdf_1y)<-c("Country","Year","Any Natural Disaster",
                       "Country Year Fixed Effects","Prob of Conflict",
                       "Predicted Conflict","Actual Conflict")
LYmatchdf_1yr<-na.omit(data.frame(compdf_1y[compdf_1y$`Predicted Conflict`==1 & compdf_1y$`Actual Conflict`==1,]))
LYmatches_1yr<-nrow(LYmatchdf_1yr)
#False neg/false positive
falsenegdf_1yr<-na.omit(data.frame(compdf_1y[compdf_1y$`Predicted Conflict`==0 & compdf_1y$`Actual Conflict`==1,]))
falseneg_1yr<-nrow(falsenegdf_1yr)
falseposdf_1yr<-na.omit(data.frame(compdf_1y[compdf_1y$`Predicted Conflict`==1 & compdf_1y$`Actual Conflict`==0,]))
falsepos_1yr<-nrow(falseposdf_1yr)
overall_1yr<-data.frame(LYmatches_1yr,falseneg_1yr,falsepos_1yr)
overall_1yr[nrow(overall_1yr) + 1,] <- c(round(LYmatches_1yr/totwars*100,2),round(falseneg_1yr/totwars*100,2),round(falsepos_1yr/totwars*100,2))
colnames(overall_1yr)<-c("Matches","False Negatives","False Positives")
rownames(overall_1yr)<-c("Number","Percentage")

pred_df_3yr$War = ifelse(pred_df_3yr$`Prob of Conflict` > .172, 1, 0)
compdf_3y<-data.frame(pred_df_3yr,totevents$War)
colnames(compdf_3y)<-c("Country","Year","Any Natural Disaster",
                       "Country Year Fixed Effects","Prob of Conflict",
                       "Predicted Conflict","Actual Conflict")
LYmatchdf_3yr<-na.omit(data.frame(compdf_3y[compdf_3y$`Predicted Conflict`==1 & compdf_3y$`Actual Conflict`==1,]))
LYmatches_3yr<-nrow(LYmatchdf_3yr)
#False neg/false positive
falsenegdf_3yr<-na.omit(data.frame(compdf_3y[compdf_3y$`Predicted Conflict`==0 & compdf_3y$`Actual Conflict`==1,]))
falseneg_3yr<-nrow(falsenegdf_3yr)
falseposdf_3yr<-na.omit(data.frame(compdf_3y[compdf_3y$`Predicted Conflict`==1 & compdf_3y$`Actual Conflict`==0,]))
falsepos_3yr<-nrow(falseposdf_3yr)
overall_3yr<-data.frame(LYmatches_3yr,falseneg_3yr,falsepos_3yr)
overall_3yr[nrow(overall_3yr) + 1,] <- c(round(LYmatches_3yr/totwars*100,2),round(falseneg_3yr/totwars*100,2),round(falsepos_3yr/totwars*100,2))
colnames(overall_3yr)<-c("Matches","False Negatives","False Positives")
rownames(overall_3yr)<-c("Number","Percentage")

pred_df_5yr$War = ifelse(pred_df_5yr$`Prob of Conflict` > .172, 1, 0)
compdf_5y<-data.frame(pred_df_5yr,totevents$War)
colnames(compdf_5y)<-c("Country","Year","Any Natural Disaster",
                         "Country Year Fixed Effects","Prob of Conflict",
                         "Predicted Conflict","Actual Conflict")
LYmatchdf_5yr<-na.omit(data.frame(compdf_5y[compdf_5y$`Predicted Conflict`==1 & compdf_5y$`Actual Conflict`==1,]))
LYmatches_5yr<-nrow(LYmatchdf_5yr)

falsenegdf_5yr<-na.omit(data.frame(compdf_5y[compdf_5y$`Predicted Conflict`==0 & compdf_5y$`Actual Conflict`==1,]))
falseneg_5yr<-nrow(falsenegdf_5yr)
falseposdf_5yr<-na.omit(data.frame(compdf_5y[compdf_5y$`Predicted Conflict`==1 & compdf_5y$`Actual Conflict`==0,]))
falsepos_5yr<-nrow(falseposdf_5yr)
overall_5yr<-data.frame(LYmatches_5yr,falseneg_5yr,falsepos_5yr)
overall_5yr[nrow(overall_5yr) + 1,] <- c(round(LYmatches_5yr/totwars*100,2),round(falseneg_5yr/totwars*100,2),round(falsepos_5yr/totwars*100,2))
colnames(overall_5yr)<-c("Matches","False Negatives","False Positives")
rownames(overall_5yr)<-c("Number","Percentage")

print(xtable(overall_1yr, type = "latex"), file = "/Users/roniyadlin/Dropbox/Climate conflict/output/tables/LYtable_1yr.tex")
print(xtable(overall_3yr, type = "latex"), file = "/Users/roniyadlin/Dropbox/Climate conflict/output/tables/LYtable_3yr.tex")
print(xtable(overall_5yr, type = "latex"), file = "/Users/roniyadlin/Dropbox/Climate conflict/output/tables/LYtable_5yr.tex")

##Comparison with N/R
#Input N/R 2010-2020 data
NR_pred<-read.csv("/Users/roniyadlin/Dropbox/Climate conflict/data/NRdata.csv")
NRroc<-roc(NR_pred$Actual.War,NR_pred$Civil_War_prob)
#N/R war prediction
NR_pred$War = ifelse(NR_pred$Civil_War_prob > .172, 1, 0)
NRtot<-sum(na.omit(NR_pred$Actual.War))
NRmatchdf<-na.omit(data.frame(NR_pred[NR_pred$War==1 & NR_pred$Actual.War==1,]))
NRmatches<-nrow(NRmatchdf)
falsenegdf_NR<-na.omit(data.frame(NR_pred[NR_pred$War==0 & NR_pred$Actual.War==1,]))
falseneg_NR<-nrow(falsenegdf_NR)
falseposdf_NR<-na.omit(data.frame(NR_pred[NR_pred$War==1 & NR_pred$Actual.War==0,]))
falsepos_NR<-nrow(falseposdf_NR)
overall_NR<-data.frame(NRmatches,falseneg_NR,falsepos_NR)
overall_NR[nrow(overall_NR) + 1,] <- c(round(NRmatches/NRtot*100,2),round(falseneg_NR/NRtot*100,2),round(falsepos_NR/NRtot*100,2))
colnames(overall_NR)<-c("Matches","False Negatives","False Positives")
rownames(overall_NR)<-c("Number","Percentage")

#Comparison between 3 ROC curves
gl_LY_comp<-ggroc(list(Year1=roc1, Year3=roc3, Year5=roc5),legacy.axes = TRUE)
AUCcomp<-data.frame(roc1$auc,roc3$auc,roc5$auc)
glsave_comp<-gl_LY_comp + xlab("FPR") + ylab("TPR") +
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="blue", linetype="dashed")
ggsave("/Users/roniyadlin/Dropbox/Climate conflict/output/figures/ROC_comp_years.png",glsave_comp,device="png")

gl_comp_all<-ggroc(list(LY_1Year=roc1, LY_3Year=roc3, LY_5Year=roc5, NR=NRroc),legacy.axes = TRUE)
AUCcompall<-data.frame(roc1$auc,roc3$auc,roc5$auc,NRroc$auc)
glsave_comp_all<-gl_comp_all + xlab("FPR") + ylab("TPR") +
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), color="blue", linetype="dashed")
ggsave("/Users/roniyadlin/Dropbox/Climate conflict/output/figures/ROC_comp_all.png",glsave_comp_all,device="png")

#Map
#Data for LY Map
mapLY5yr<-data.frame(compdf_5y,predictdf$latitude,predictdf$longitude)
mappredlat<-na.omit(data.frame(mapLY5yr$predictdf.latitude[mapLY5yr$Predicted.Conflict==1]))
mappredlong<-na.omit(data.frame(mapLY5yr$predictdf.longitude[mapLY5yr$Predicted.Conflict==1]))
mappred<-data.frame(mappredlat,mappredlong)
colnames(mappred)<-c("Lat","Long")
write.csv(mappred,"/Users/roniyadlin/Dropbox/Climate conflict/output/figures/Mapping/predictLY.csv",row.names = FALSE)
mapactlat<-na.omit(data.frame(mapLY5yr$predictdf.latitude[mapLY5yr$Actual.Conflict==1]))
mapactlong<-na.omit(data.frame(mapLY5yr$predictdf.longitude[mapLY5yr$Actual.Conflict==1]))
mapact<-data.frame(mapactlat,mapactlong)
colnames(mapact)<-c("Lat","Long")
write.csv(mapact,"/Users/roniyadlin/Dropbox/Climate conflict/output/figures/Mapping/actualLY.csv",row.names = FALSE)

#Data for Matches

matchdf<-na.omit(data.frame(mapLY5yr[mapLY5yr$Predicted.Conflict==1 & mapLY5yr$Actual.Conflict==1,]))
mapmatch<-data.frame(matchdf$predictdf.latitude,matchdf$predictdf.longitude)
colnames(mapmatch)<-c("Lat","Long")
write.csv(mapmatch,"/Users/roniyadlin/Dropbox/Climate conflict/output/figures/Mapping/matchLY.csv",row.names = FALSE)

#mapmatch<-data.frame(LYmatchdf_5yr,predictdf$latitude,predictdf$longitude)
#mappredlat<-na.omit(data.frame(mapLY5yr$predictdf.latitude[mapLY5yr$Predicted.Conflict==1]))
#mappredlong<-na.omit(data.frame(mapLY5yr$predictdf.longitude[mapLY5yr$Predicted.Conflict==1]))
#mappred<-data.frame(mappredlat,mappredlong)
#colnames(mappred)<-c("Lat","Long")
#write.csv(mappred,"/Users/roniyadlin/Dropbox/Climate conflict/output/figures/Mapping/predictLY.csv",row.names = FALSE)

#NRmappredlat<-na.omit(data.frame(NR_pred$Latitude[NR_pred$War==1]))
#NRmappredlong<-na.omit(data.frame(NR_pred$Longitude[NR_pred$War==1]))
#NRmappred<-data.frame(NRmappredlat,NRmappredlong)
#colnames(mappred)<-c("Lat","Long")
#write.csv(mappred,"/Users/roniyadlin/Dropbox/Climate conflict/output/figures/Mapping/predictNR.csv",row.names = FALSE)

